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ABSTRACT 

Astrophysical observations originate from matter that interacts with radiation or transported parti- 
cles. We develop a pragmatic approximation in order to enable multi-dimensional simulations with 
basic spectral radiative transfer when the available computational resources are not sufficient to solve 
the complete Boltzmann transport equation. The distribution function of the transported particles 
is decomposed into a trapped particle component and a streaming particle component. Their sepa- 
rate evolution equations are coupled by a source term that converts trapped particles into streaming 
particles. We determine this source term by requiring the correct diffusion limit for the evolution of 
the trapped particle component. For a smooth transition to the free streaming regime, this 'diffu- 
sion source' is limited by the matter emissivity. The resulting streaming particle emission rates are 
integrated over space to obtain the streaming particle flux. Finally, a geometric estimate of the flux 
factor is used to convert the particle flux to the streaming particle density, which enters the evalua- 
tion of streaming particle-matter interactions. The efficiency of the scheme results from the freedom 
to use different approximations for each particle component. In supernovae, for example, reactions 
with trapped particles on fast time scales establish equilibria that reduce the number of primitive 
variables required to evolve the trapped particle component. On the other hand, a stationary-state 
approximation considerably facilitates the treatment of the streaming particle component. Differ- 
ent approximations may apply in applications to stellar atmospheres, star formation, or cosmological 
radiative transfer. We compare the isotropic diffusion source approximation with Boltzmann neu- 
trino transport of electron flavour neutrinos in spherically symmetric supernova models and find good 
agreement. An extension of the scheme to the multi-dimensional case is also discussed. 
Subject headings: supernovae: general — neutrinos — radiative transfer — hydrodynamics- -methods: 
numerical 



1. INTRODUCTION 

Most applications in computational astrophysics in- 
volve matter that is in thermal or reactive equilibrium. 
Astrophysical observations, however, require information 
to propagate from the astrophysical site to a terrestrial 
observer. Hence, the observationally interesting events 
involve additional particle species that are not trapped 
in the fluid and must therefore be treated by radiative 
transfer. In many traditional numerical models, the equi- 
librated matter is treated in the hydrodynamic limit, 
while transported particle species are treated by a suit- 
able algorithm of radiative transfer. If this splitting in 
the numerical method is based on the particle species, it 
can become extremely challenging and inefficient if the 
range of conditions in the astrophysical scenario is large. 
In situations where radiative particles are in thermal or 
reactive equilibrium, the non-local algorithm of radiative 
transfer has to be perfectly consistent with the hydrody- 
namics scheme and capable of handling very stiff source 
terms in order to evolve the applicable equilibria. 

For example core collapse supernova events 
have cha llenged compute r mod els for s everal 
decades (iColgate fc White! U963 lArnettl \\M& 



Bowers fc Wilson! Il982t iBethe fc Wilson! Il985t fBruenrJ 
19851: iMvra fc Bludmar]ll989h . These stellar explosions 
occur at the end of the life of a massive star when the 
growing iron core, as the end product of nuclear fusion, 
becomes unstable against gravitational collapse. The 
collapse is halted only after nuclear density is reached 
at the centre. The repulsive strong interaction and the 



neutron degeneracy pressure reduce the compressibility 
of matter so that a sound wave travels outward through 
the causally connected inner core. It turns into a shock 
wave when it reaches supersonically accreting matter 
from the outer layers and, due to dissociation and 
neutrino losses, stalls within few milliseconds to an 
initially hydrostatic accretion front, which is thought 
to slowly expand over several hundreds of milliseconds 
before the explosion sets in. The new-born protoneutron 
star (PNS) at the centre of the event has a larger density 
by many orders of magnitude than the shock-heated 
and dissociated matter that continues to accumulate 
behind the accretion front. This accreted matter is 
subject to a variety of fluid instabilities and asymmetric 
flow patterns that couple to the energetic neutrino 
radiation field produced within or close to the surface 
of the PNS. The emergence of an explosion and the 
detailed interaction between the multi-dimensional flow 
and the neutrino field is subject to a long-standi ng and 
ongoing debate. See for example (iBethd I1990D for a 
theor e tical description of t h e scenario or (iBruenn et al.l 
[20061 iDessart et all 120071: iMarek fc Jankal 120071) for 
recent axisymmetric computer simulations with spectral 
neutrino transport. 

Several aspects of the core collapse supernova scenario 
make computer models difficult to perform: The rele- 
vant densities in the computational domain range over 
ten decades from about 10 15 g/cm 3 at the centre of the 
PNS down to about 10 5 g/cm 3 in the outermost lay- 
ers (the star is much larger but does not dynamically 
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react to the changes at the centre during the first few 
hundred milliseconds after bounce). This density range 
implies very disparate dynamical time scales. At high 
densities, the weak interaction rates are much faster than 
the dynamical time scale, while at low density, they are 
even slower than the evolution time scale. This requires 
implicit finite differencing of many processes to perform 
simulations in a reasonable amount of time. Moreover, 
the gravitational and stored internal energy scales are of 
order 10 53 erg, while typical kinetic explosion energies 
are two orders of magnitude smaller. The presence of 
a strong shock with near-relativistic infall velocities fur- 
ther challenges the resolution and stability of numerical 
models. In this very dynamical environment with multi- 
dimensional fluid instabilities and differential rotation, a 
computer model has to quantify the exchange of energy 
between matter and the emitted neutrino field in order 
to allow statements about the explosion dynamics to be 
made. As the neutrino cross sections at given matter 
conditions scale with the square of the neutrino energy, 
the neutrino transport problem must be considered as a 
superposition of several single-energy (monochromatic) 
transport problems, where the geometry of the scattering 
sphere depends on the neutrino energy and type. This 
leads to multi-group transport schemes where all propa- 
.ation quantitie s are solved i n separate energy groups 
Mazurekl 119751 : lArnettl Il977| ). Moreover, stationary- 
state studies with spherically symmetric Boltzmann neu- 
trino transport showed that the angular distribution of 
the neutrino propagation directions have also to be taken 
into acount to obtain accurate neutrino heating and 
cooling rates |WilsorJll97lUMezzacappa fe Bruennll 19931: 
Messer et al.l Il998l lYamada et all 119991 : iBurrows et al.l 
2003). 



In the mean time, the spherically symmetric Boltz- 
mann transport equation is routinely solved for all 
neutrino flavours in dynamical mo dels of stellar core 
collap s e and postbounc e evolu tion (iLiebendorfer et al.l 
20011: iRampp fe Jankal 120021 : iThompson et all 120031 : 



Sumivoshi et al.ll2005l ) . These comprehensive and sophis- 
ticated codes treat disparate time scales by implicit finite 
differencing and require large computational resources. 
However, due to the enormous density range in the sce- 
nario, a considerable fraction of the computational do- 
main is either neutrino-opaque or neutrino-transparent. 
In a fully neutrino-opaque regime, a numerical solution 
to the Boltzmann transport equation is rarely more ac- 
curate than the solution of the simpler diffusion equation 
that describes the correct physical limit. Similarly, the 
angular distribution of the neutrino propagation direc- 
tions in the far-field free streaming regime is precisely 
determined by the geometry and emissivity of the dis- 
tant sources, while an angular discretisation of the Boltz- 
mann equation may lead to discretisation errors. Hence, 
we suggest that a combination of different approaches 
to radiative transfer is applied in the same astrophysi- 
cal simulation. This leads to the concept of an 'adaptive 
algorithm', which should not be regarded as an incon- 
sistency or drawback. On the contrary, it saves compu- 
tation power that can be invested for better resolution, 
more input physics, or the coverage of a larger param- 
eter space. Current state-of-the-art supernova models 
with spectral neutrino transport are performed under 
the assumption of axisymmetry. As a 2D solution of 



the Boltzmann equation is computationally inefficient for 
the reasons discussed above, different groups either ap- 
ply ID Boltzmann solutions in separate angular segments 
in combination with 2D neutrino advection in opaque 
regimes (ray- by-ray) (|Marek fc Jankal l2007f ). resort to 
the multi -group flux-limited diffusion (MGFL D) approx- 
imation (|Walder et all 120051: lOtt et al.l I2008T) . or use a 
combination of both (|Bruenn et al.ll2006i r 

Three-dimensional supernova models were hith- 
erto only affordable with 'grey' neutrino transport 
(|Frver fe W arren 2004) or with other simple ad hoc ap- 
proximations to the neutrino physics (IScheck et al.ll2004l : 
lOtt et alJ 120071 IScheidegger et al.ll2C)08t ). But many fea- 
tures of the supernova should be studied in three spa- 
tial d imensions, for exam ple the pattern of the accretion 
flow dHerant et al.lll994h. the s tanding accretion shock 
instability (Blondin et al. 2003), the excitation of PNS 
g-modes (|Burrows et al.ll2006D . and the ev olution of mag- 
netic fields (see e.g. iKotake et aTl (|2006f ) and references 
therein). In order to enable three-dimensional supernova 
models with spectral neutrino transport and in order to 
accelerate parameter studies with simulations of lower di- 
mensionality, we construct the isotropic diffusion source 
approximation based on several ideas scattered across 
the literature. Our algorithm is not meant to compete 
in accuracy with more detailed solutions of the transport 
equation in the semi-transparent regime. As the transi- 
tion from the opaque to the transparent regime will be 
handled by interpolation and basic geometric considera- 
tions, it is important that the approach is verified by an 
accurate transport solution for every new field of applica- 
tion. Our approach shares this limitation with the well- 
known flux-limited diffusion approximation. The main 
advantage of the isotropic diffusion source approxima- 
tion over the latter is that the fluxes and flux factors 
in the transparent regime are determined by the non- 
local distribution of sources rather than the local inten- 
sity gradient, which may give an incorrect flux direction. 
Moreover, the new approach tries to limit the computa- 
tionally expensive solution of the multi-dimensional dif- 
fusion equation to the opaque regime where the diffu- 
sion approximation is adequate. The goal is to create 
a flexible algorithm that efficiently implements only the 
dominant features of radiative transfer in one consistent 
framework. 

The most obvious method of avoiding the inefficient 
global application of algorithms is a decomposition of the 
problem in space, so that one algorithm is used in one 
subdomain while another algorit hm is used in another 
subdomain (jChick fe Cassenl Il997h . Another approach 
is the decomposition in the momentum phase space into 
particles with thermal velocities and particles with supra- 
thermal velocities. If one treats the thermal particles by 
an efficie nt fluid model, a hybri d kinetic/fluid model is 
obtained (|Crouseilles et all 2004) . The isotropic diffusion 
source approximation described in this article is more 
similar to the so-cal led <5/-method (Par ker fe Led [19931 : 
iBrunner et al.l fl999h . We also decompose the distribu- 
tion function of the transported particles into a ther- 
mal component and a perturbation that is allowed to 
overlap the thermal part in the entire computational do- 
main and particle phase space. However, in our case, 
the perturbation is not considered to be small. Inspire d 
by flux-limited diffusion (|Levermore fe PomraningTl 981). 
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the evolution of the radiation component is guided by the 
diffusion limit at large opacities and by the free streaming 
limit at low opacities. In contrast to flux- limited diffu- 
sion, we build our scheme conceptually on approxima- 
tions of the collision integral rather than particle fluxes. 
We assume that the free streaming particle flux is dom- 
inated by the flux emerging from the diffusive domain 
so that the sources for the far field can easily be inte- 
grated. In the stationary-state limit, the determination 
of the par ticle fluxes reduces to the solution of a Poisson 
equation (|Gnedin fe Abe]||2001h . Since the matter inter- 
acts according to the local particle abundances and not 
fluxes, we have to convert the particle fluxes to local par- 
ticle densities. This is achieved by a geometric estimate 
of the flux factor as suggest ed and evaluated by Bruenn 
in ([Liebendorfer et al.ll2004f ). 

In Sect. [5] we describe in detail how these concepts 
enter the framework of the isotropic diffusion source ap- 
proximation, which we design for the transport of mass- 
less fermions through a compressible gas. Its connection 
to the well-known diffusion limit is made in Appendix A. 
In Sect. [3l we evaluate the performance of this approxi- 
mation in comparison with Boltzmann neutrino trans- 
port in spherical symmetry. Finally, in Sect. |U we 
discuss the extension to multi-dimensional simulations. 
Details of the finite differencing and implementation are 
given in Appendix B. 

2. THE ISOTROPIC DIFFUSION SOURCE 
APPROXIMATION (IDSA) 

In the IDSA, the separation into hydrodynamics and 
radiative transfer is not based on particle species, but on 
the local opacity. One particle species is allowed to have a 
component that evolves in the hydrodynamic limit, while 
another component of the same particle species is treated 
by radiative transfer. The restriction of a chosen radia- 
tive transfer algorithm to the more transparent regimes 
enables the use of more efficient techniques that would 
not be stable in the full domain. In opaque regimes, on 
the other hand, one can take advantage of equilibrium 
conditions to reduce the number of primitive variables 
that need to be evolved. This algorithmic flexibility can 
drastically decrease the overall computational cost with 
respect to a traditional approach. 

In the IDSA, we decompose the distribution function 
of one particle species, /, into an isotropic distribution 
function of trapped particles, and a distribution func- 
tion of streaming particles, / s . In terms of a linear op- 
erator D () describing particle propagation, the trans- 
port equation is written as D (/ = /' + / s ) = C, where 
C = C* + C s is a suitable decomposition of the collision 
integral according to the coupling to the trapped (C*) or 
streaming (C s ) particle components. The ansatz 



D (/ t )=C t - £ 
£(f) = C s + £, 



(1) 
(2) 



requires that we specify an additional source term S, 
which converts trapped particles into streaming parti- 
cles and vice versa. We determine it approximately from 
the requirement that the temporal change of / in Eq. 
(fTj) has to reproduce the diffusion limit in the limit of 
small mean free paths. Hence, we call £ the 'diffusion 
source'. In regions of large mean free paths, we limit the 
diffusion source by the local particle emissivity. Once £ 



is determined by the solution of Eq. (fTJ) for the trapped 
particle component, we calculate the streaming particle 
flux according to Eq. © by integrating its source, C s +£, 
over space. Finally, the streaming particle distribution 
function / s is determined from the quotient of the net 
particle flux and a geometric estimate of the flux factor. 
The diffusion source will turn out to have an additional 
weak dependence on / s . Thus, iterations or information 
from past time steps will be used in the above sequence 
to reach a consistent solution. 

2.1. Application to radiative transfer of massless 
particles 

As our target application is neutrino transport in 
core collapse supernovae, we develop and test the IDSA 
using the example of th e Q(v/c) Boltzmann equation 
in spherical symmetry dLindauisd Il966t ICastorl 119721 : 
Mezzaca ppa fe Bruennlll993lK 
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This transport equation describes the propagation of 
massless fermions at the speed of light, c, with respect 
to a compressible background matter having a rest mass 
density p. The particle distribution function / (t, r, p, E) 
depends on the time, t, radius, r, and the momentum 
phase space spanned by the angle cosine, p, of the parti- 
cle propagation direction with respect to the radius and 
the particle energy, E. The momentum phase space vari- 
ables are measured in the frame comoving with the back- 
ground matter, which moves with velocity v with respect 
to the laboratory frame. We denote the Lagrangian time 
derivative in the comoving frame by df/dt. Note that 
the derivatives df /dp and df /dE in Eq. ([3]) are also 
understood to be taken comoving with a fluid clement. 
The particle density is given by an integration of the 
distribution function over the momentum phase space, 
n(t,r) = 47t/ (hef J f (t, r, p, E) E 2 dEdp, where h de- 
notes Plancks constant. On the right hand side, we in- 
clude a particle emissivity, j, and a particle absorptiv- 
ity, x, as well as an isoenergetic scattering kernel, R. 
We write out all blocking factors (1 — /) in Eq. ([3]) to 
ease the identification of in-scattering and out-scattering 
terms. The shorthand notation /' refers to / (t, r, p! ,E), 
where p' is the angle cosine over which the integration is 
performed. For the present state of our approximation 
we neglect inelastic scattering. 

2.2. Trapped particles 

We separate the particles described by the distribution 
function / = /* + / s in Eq. ^ into a 'trapped particle' 
component, described by a distribution function / , and 
a 'streaming particle' component, described by a distri- 
bution function / s . We assume that the two components 
evolve separately according to Eq. ([3]), coupled only by 
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an as yet unspecified source function E which converts 
trapped particles into streaming ones and vice versa. In 
this subsection we discuss the evolution equation of the 
trapped particle component, 
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We assume that the distribution of the trapped particle 
component, /* = /* (t, r, E), and the source function, E, 
are isotropic. The angular integration of Eq. |4| then 
reduces to 



df Idhxp df 
cdt 3 cdt dE 



: 3 - (3 +X)f~ S. (5) 



However, even if we are now steering towards the hy- 
drodynamic limit, the evolution of the trapped parti- 
cle distribution function in our approximation should at 
least reproduce the correct diffusion limit. The physical 
understanding of diffusion relies on fast-moving parti- 
cles with a very short transport mean free path (see the 
derivation of the diffusion limit, Eq. (|39p . in Appendix 
A). The divergence of the small net particle flux leads 
to a slow drain (or replenishment) of particles in a fluid 
element. In order to accomodate this diffusive drain (or 
replenishment) of trapped particles in Eq. ([5]), we have to 
implement it through the so far unspecified source term 
E. In the framework of our approximation, trapped par- 
ticles are converted to streaming ones (or vice versa) at 
the same rate diffusion would drain (or replenish) parti- 
cles in the fluid element. A comparison of Eq. ([S]) with 
Eq. ([39]) in Appendix A suggests the following diffusion 
source, 



10/ - r 2 df 
r 2 dr \ 3 (j + x + 4>) dr 



+ (3 + x) \ I r<h'- 



(6) 

Isoenergetic scattering enters Eq. ([5]) only via its opac- 
ity <j) in the expression for the transport mean free path 
A = 1/ (j + \ + <t>) that determines the flux in Eq. ([6]) 
(see Appendix A). The additional term {j + %) /2j fdp 
accounts for the absorption of streaming particles in mat- 
ter. Its necessity is best understood with the help of Fig. 
([T]). Shown are the different particle fluxes that relate 
to one fluid element. The fluid element contains both 
matter and trapped particles. They interact by the emis- 
sivity j and absorption (j + x) f (vertical arrows). The 
radiation particles leaving the fluid element will convert 
into streaming particles at the rate E (upper horizon- 
tal arrow). We do not include the possibility of a direct 
emission of matter into the streaming particle compo- 
nent because this process is more easily accounted for 
by a high conversion rate E of trapped particles. On 
the other hand, we account for streaming particles that 
are absorbed by matter as represented by the arrow de- 
noted by (j + x) f s m Fig- ©■ ln the diffusion limit, 
the net particle exchange of the fluid element with its 
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Fig. 1. — The shaded box schematically represents a fluid element 
in the diffusion source approximation. It contains matter (lower 
part) and trapped radiation particles (upper part). The interaction 
with other fluid elements can exclusively occur by the exchange of 
streaming particles or the combined hydrodynamics of the matter 
and trapped particles. Hence, streaming particles can be absorbed 
in matter at the rate (j + \) f B and trapped particles are converted 
to streaming particles at the rate S. Within the fluid element, 
matter emits trapped particles at the rate j and absorbs trapped 
particles at the rate (j' + x)/*- The emissivity in the absorption 
term originates from the identity j (1 — /) — x/ = 3 ~~ (j + x) /> 
which hides the Pauli blocking factor in the absorption term. 



f s dp, must correspond to 
I). This is the case for our 



environment, E — (J + x) /2 
the diffusion term in Eq. ([3! 
choice of E in Eq. ©. 

If the above scheme is applied in a more transparent 
regime where the diffusion approximation does not hold, 
the diffusion source in Eq. ((6|) may become arbitrarily 
large. This would be inconsistent with the particle fluxes 
drawn in Fig. JT]) because over a long time the diffusion 
source can not exceed the emissivity of trapped particles 
without creating an unphysical deficit in trapped parti- 
cles. Instead, we would like the trapped particle compo- 
nent to drop to zero and stay zero in the limit of long 
mean free paths. This is achieved if we limit the diffu- 
sion source in Eq. © to E < j. If the diffusion source 
and emissivity reach equality, the matter absorptivity - 
(j + x) f removes remaining trapped particles while all 
newly emitted ones are directly converted to streaming 
particles that escape the fluid element. With this limit 
imposed, the net interaction of particles with matter in 
Fig. ([I]) has the correct limit for large mean free paths, 

3 - (3 + X) f- 

The search for a lower bound to the diffusion source 
is less straight-forward. In principle, a negative E can 
not be physically excluded. It corresponds to stream- 
ing particles that become trapped in a region of large 
opacity and low absorptivity. A limit corresponding to 
the restriction f < 1 would suggest E > — x- How- 
ever, / = 1 can be an excessively large particle density 
compared to the physically expected value. Much more 
stable and more accurate results were obtained by the 
requirement that /* < j / (j + x), where j / (j + x) rep- 
resents the equilibrium distribution function. If one con- 
siders this condition for the particle fluxes in Fig. {T]) 
this leads to the simple requirement E > 0. Hence, the 
net absorption of particles in a fluid element can not ex- 
ceed (j + x) f s ■ The diffusion source from Eq. ([6]) then 
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becomes 



mm < max 



a + (j + X) 

J2, 



1 

2 



fdp,0 



1 d 



r 2 dr \ 3 ( j ' + x + 4 1 ) dr 



(7) 



Assuming that the streaming particle density 1/2 J f s dp 
was known from the considerations described in the fol- 
lowing subsection, Eq. can be used to calculate the 
evolution of the trapped particle component consistently 
with the diffusion source specified in Eq. ([7]). 

2.3. Streaming particles 

The evolution equation for the streaming particle com- 
ponent consists of all terms in Eq. ([3]) that have not been 
considered in the evolution equation ^ for the trapped 
particles. As mentioned above, we neglect direct scat- 
tering from the trapped component into the streaming 
component and vice versa. That is 
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For the evolution of streaming particles we neglect the 
scattering integrals on the right hand side of Eq. ©, 
because the streaming particle density is designed to be 
small compared to the trapped particle density in re- 
gions where the scattering rate dominates the interac- 
tions. The weak coupling between streaming particles 
and the background matter in more transparent regimes 
where the streaming particle component becomes large 
makes an inertial laboratory frame more convenient for 
the solution of Eq. ((8|) than the frame comoving with 
the fluid. One may simply substitute the conditions 
din pi cdt = and v = for a static background in Eq. 
|8]l to find the streaming particle evolution equation in 
the laboratory frame: 

Note that the quantities carrying a hat now are mea- 
sured in the laboratory frame. In this frame, the particle 
energy is a constant of motion. The cumbersome energy 
derivatives of the distribution function in Eq. (jBJ) are 
avoided, but some quantities require Lorentz transfor- 
mations between the comoving and laboratory frames. 

In contrast to the full Boltzmann equation (J3j> , the 
source term on the right hand side couples only weakly 
to the particle distribution function / s . In the diffusive 
domain, the angle integrated term — (j + j^J f s cancels 
to a large extent with the corresponding contribution to 
£ as defined in Eq. ([7]). In the free streaming domain, 

the emissivity j and opacity x are assumed to be small. 
Depending on the application, there might be various 
suitable approaches to solve Eq. ©. Here, we will pro- 
ceed with a stationary-state approximation. If we assume 



that the fluxes of free streaming particles reach station- 
ary values much faster than the dynamical time scale 
of interest, we can drop the time-derivative in the first 
term in Eq. ([9]). If the source on the right hand side is 
assumed to be known from a consistent solution of Eqs. 
((5|) and jJJ, the integration of Eq. (J9j) over angles leads 
to a Poisson equation for a potential ip, whose gradient 
represents the particle flux: 

d * 1 tfm 



_L__9 

r 2 dr 



dr 

t chp\ 
dr J 



- i+% / s +e dp. (io) 



In the end, we are interested in the streaming parti- 
cle density rather than flux. Both quantities are related 
by the flux factor. It is a very challenging problem to 
calculate an accurate flu x factor. However, Bruenn in 
(jLiebendorfer et all l2004T l suggested a simple and very 
successful approximation based on the assumption that 
all particles of a given energy are isotropically emitted at 
their corresponding scattering sphere. This leads to the 
additional relation 



~ / nEuip 



J dr 



1 



(11) 



1 



max(r\ /?.„(£)) 



where R v (E) is the radius of the monochromatic scat- 
tering sphere that depends on the particle energy E. 

In order to solve for both trapped and streaming par- 
ticle components, we have to Lorentz-transform quanti- 
ties between the comoving and laboratory frames. The 
source of Eq. (fTOf in the laboratory frame is related to 
the comoving frame by 

' 1 \-(j + x)f s + ^ 



dpE z dE 



[- (j + X )f s + n dpE 2 dE. 



(12) 



Here we used the fact that the source is isotropic 
in the comoving frame together with the relations 
dpE dE = dpEdE, E = -f(l + pv/ c)E and jE = 
jE (|Mihalas & Weibel Mihalasl Il984f) . where 7 = 

1/ *Jl — (u/c) 2 is the Lorentz factor. On the other hand, 
the streaming particle contribution to Eq. J7]) in the co- 
moving frame is related to the solution of Eq. (jTUJ) in the 
laboratory frame by 

i f fdpE 2 dE = X -l /» 7 (1 - pv/c) dpE 2 dE 



= 7 



(13) 



where the first term is given in Eq. pip . In O (v/c) we 
may neglect the factor 7. 

Hence, as soon as the evolution of the trapped parti- 
cle component according to Eq. (O has determined the 
source in Eq. (7|), one can transform it to the laboratory 
frame by Eq. (fT2"|) . Then one determines the stationary- 
state streaming particle flux and density based on Eqs. 
(jTUJ) and (fTTj) in the laboratory frame. The resulting 
streaming particle density is transformed back to the co- 
moving frame by Eq. (|13p . In this form it can be used 
for the next iteration or time step in Eq. ([7]). 
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2.4. Coupling with hydrodynamics 

The evolution of the trapped particle component in 
Eqs. ([5]) and ([7]) and the stationary-state limit of the 
streaming particle component in Eqs. (jlOp and (llip must 
be coupled with the dynamics of the background matter. 
We assume that the dynamics of the background matter 
is well described by the conservation laws of hydrody- 
namics, 



0_ 

di 



r z dr v ' 



(14) 



where U is a vector of primitive variables and F a vec- 
tor of fluxes. We rewrite the evolution equation ([5]) for 
the trapped particle component with an Eulerian time 
derivative and use the continuity equation to substitute 
the dhip/dt term by the velocity divergence, 
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(15) 



One could now merge the hydrodynamical update on the 
left hand side of Eq. (fl"5|) with the conservation law in 
Eq. (fT4"| . However, if the dynamics of the background 
matter requires a large-scale numerical simulation and 
if there are many particle species, the required memory 
and CPU time to store and advect all these quantities can 
easily exceed the capacity of the available hardware. In 
our supernova application, for example, we need in three 
dimensions a vector U of 6 entries to describe the hy- 
drodynamics as detailed below. The neutrino transport 
involves 2-4 neutrino species (v c , D c , ^ t / r , ^Vt) with at 
least 12 entries per species to sample the different neu- 
trino energies E in /* (t, r, E). 

If the characteristic time scale of local reactions be- 
tween the transported particle species is faster than the 
diffusion time scale, one can use equilibrium conditions 
to reduce the number of primitive variables necessary to 
describe the distribution functions of trapped particles. 
In the supernova model, for example, we can approxi- 
mate the spectrum by a thermal spectrum. Note, that 
this assumption is only made for the trapped particles 
within a fluid element. The streaming particles, which 
communicate between fluid elements, keep their detailed 
spectral information. Our scheme should therefore not 
be confused with a 'grey' scheme that may additionally 
assume a predefined energy-spectrum for the particle ex- 
change between fluid elements. We characterise the ther- 
mal equilibrium by a particle number fraction, Y*, and 
a particle mean specific energy, Z , 



?7l b 47T 

P (hef 
m\, 47r 



fE 2 dEdp 
f t E 3 dEdfi, 



(16) 



where m b is a constant relating a particle number to the 
rest mass of the background matter (in our application 
it is the baryon mass). 

The corresponding energy integrals performed on Eq. 
(fT5|) lead to evolution equations for y* and Z t ) 



d_ 

di 



(pY* 







r2 dr 



(^vpY 1 ) 



-m h - 



4nc 
(hef 



J - (J + X) f 
8 



T,dp 



E 2 dE 
pZ* 



= m b 



4-7TC 

{hef 



I [j - (j + X)f- E] E 3 dE. 



The equation for Z t corresponds to an energy equation 
with a pdV term for the radiation pressure pZ t /3 on the 
left hand side. Analogously to an entropy equation it 
can be cast into a conservative form for the evolution of 
(pZ*) ' . Hence, one can solve the advective part of Eq. 



15j) together with the hydrodynamics conservation law 
(|14|) based on the following primitive variables: 
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where p is the fluid pressure, e the fluid specific internal 
energy and Y c the electron fraction. The index I labels 
different species of trapped particles. 

As the distribution function of trapped particles in 
thermal equilibrium, ff (E) = {exp [0i (E — /x;)] + 1} , 
has two free parameters and pi , it can be reconstructed 
so that Eq. (fTS)) is fulfilled for the new values of and 
Z\. Then, the remaining update of the trapped particle 
distribution in Eq. (jTSj) is given by 
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This equation also determines the net interaction rates, 
si = ji — (ji + xi) (fi + ff), between matter and the ra- 
diation particles, which leads to the following changes of 
the electron fraction and internal specific energy: 
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(20) 
(21) 



The changes in the matter electron fraction and specific 
energy feed back into the emissivity and absorptivity 
used in Eq. (fl"8|) . Once a consistent solution has been 
found, Eq. (|18|) allows updates of the trapped particle 
fraction, of the trapped particle specific energy, and of 
the matter velocity, which is subject to radiation pres- 
sure. 
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Finally, the cycle of updates is completed by the solution 
of Eqs. (fT0|) and (fTT|) for the distribution function / s for 
streaming particles based on the sources determined in 
Eqs. (jT8|) and ([7]). Details on the finite differencing of 
our implementation of the IDSA are given in Appendix 
B. 

3. VERIFICATION FOR SPHERICALLY SYMMETRIC 
SUPERNOVA MODELS 

We implement the IPSA in the hydrodynamics code 
Agile (Licbendorfcr et al. 2002). For the evaluation of 
the radiative transfer approximations, we run the simu- 
lations in the Newtonian limit and include only electron 
flavour neutrinos. Their dominant emission and absorp- 
tion reactions are 

e~ + p^n + v c 

e+ + n^p + i> e , (25) 

where e~ , e + , v c , v c refer to electrons, positrons, electron 
neutrinos and electron antineutrinos, respectively, that 
interact with protons, p, and neutrons n. The opacities in 
our example are given by isoenergetic scattering on nucle- 
ons and nucl ei. All weak in teractions are implemented as 
described in (|Bruennlll985h . The thermodynamical state 
of matter as a function of density, p, temperature, T, 
and electron fraction, Y e , is calculated by the Latt imer- 
Swesty equation of state (jLattimer fc Swestvlll991h . 

Very simple and efficient approximations to treat the 
deleptonisation in the collapse phase have been suggeste d 
and evaluated in an earlier paper (jLiebendorferl 12005) . 
Here we concentrate on the postbounce phase, which is 
more difficult to capture with neutrino physics approxi- 
mations. The neutrinos not only leak out, but they also 
interact at distant locations and feed back into the hy- 
drodynamics by their transfer of lepton number and en- 
ergy. Due to the as yet restricted inp ut physics, we base 
our fi rst comparison on model Nf 3 in (Licbendorfcr ct al. 
2005), which we have repeated with the most recent 
code version to obtain perfect consistency in all parts 
that are not related to the neutrino transport. The ad- 
ditional neutrino-electron scattering and pair creation 
rates, which are included in model N13, do not signif- 
icantly contribute to the electron flavour neutrino trans- 
port in the postbounce phase. In order to keep the ap- 
proximation as simple as possible, we additionally neglect 
the transformations between comoving and laboratory 
frames in our current implementation. 

We start with the verification of the description of the 
trapped particle component. The trapped particle dis- 
tribution function is assumed to be thermal and is pa- 
rametcrised by the two parameters inverse temperature, 
(3, and the particle degeneracy parameter, fi/(kT). They 
are chosen such that Eq. (|16p is fulfilled. In Fig. [3] we 
compare these parameters to the temperature of mat- 
ter and to the equilibrium degeneracy parameter deter- 
mined by reactions (|23|) . We find as expected that the 
neutrino distribution function of a reference simulation 
with Boltzmann neutrino transport is very well described 
by the parameterisation in regimes where the neutrinos 
are trapped and thermal and weak equilibrium are good 
approximations. Deviations are visible where the neu- 
trinos start to decouple from matter. The IDSA leads 
to decreasing trapped particle abundances in this tran- 
sition regime as more and more particles change from 
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Fig. 2. — Neutrino temperature and degeneracy parameters re- 
sulting from the application of Eq. (1 16 I t at a typical time of 150 ms 
after bounce. As indicated in the legend, t he s olid lines show the 
degeneracy of the particles involved in Eq. (1 25 I t together with the 
(anti-)neutrino degeneracy parameters. The dashed lines show the 
matter temperature together with the (anti-)neutrino temperature 
parameters. 



the trapped particle component to the streaming parti- 
cle component. This is visible in the decreasing chemical 
potentials at radii larger than 40 km. The temperature 
parameter of the remaining trapped particle component 
is slightly larger than the temperature of the matter. 

This decoupling of particles from opaque matter is an 
important feature that needs to be handled by a radiative 
transfer algorithm. Since the IDSA treats this transition 
by a gradual conversion of 'trapped particles' to 'stream- 
ing particles', we can compare the sum of the trapped 
particle and streaming particle densities with the total 
particle density in a reference simulation that solves the 
Boltzmann equation. Fig. [3] shows this comparison for 
a typical time slice at 150 ms after bounce. Because 
the trapped particle distribution function in Eq. (|16l) is 
characterised by the particle abundance and mean en- 
ergy, we select these two quantities for the comparison. 
Panel (a) shows how the trapped and streaming parti- 
cle abundances overlap and add up to a reasonable total 
particle abundance. Panel (b) shows similar agreement 
in the particle mean energies. Although the results agree 
very well in general, we also note that the deviations are 
largest around 50 km radius. Fig. [^indicates that this 
is exactly in the difficult regime where the neutrinos de- 
couple from matter. The less relevant deviations at radii 
larger than 220 km stem from the Doppler energy shift 
at the shock front, that is not present in the IDSA data, 
because the transformation between the laboratory and 
comoving frames has been neglected. 

We can analyse the discrepancies further by looking at 
the particle spectra in the different regimes. The panels 
on the left hand side of Fig. 0] show neutrino spectra, 
while the panels on the right hand side show antincu- 
trino spectra. The spectra in panels (a) and (b) are 
taken at 40 km radius from the time slice at 150 ms 
after bounce. One can immediately see that the trapped 
neutrino component (dashed lines) dominates over the 
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Fig. 3. — Panel (a) shows neutrino abundances as a function of radius at a typical postbounce time of 150 ms. The thick solid line shows 
the neutrino abundance in the reference simulation with Boltzmann transport. The application of the IDSA leads to a trapped neutrino 
abundance (dashed line) and a streaming neutrino abundance (dash-dotted line). They sum up to the total neutrino abundance (thin solid 
line), which is comparable to the Boltzmann result. Graph (b) demonstrates that also the neutrino mean energies between the reference 
model (thick lines) and the IDSA (thin lines) agree nicely. 



streaming particle component (dash-dotted lines). The 
sum of both (thin solid line) matches well with the refer- 
ence solution (thick solid line). At these conditions the 
spectrum of the neutrinos is very accurately represented 
by the IDSA, while small differences become visible in the 
antineutrino spectra. The spectra match well, because 
the thermal distribution function used for the trapped 
neutrino component is still a good approximation at the 
border of the neutrino-opaque regime. Panels (c) and (d) 
show the same quantities further out, where the trapped 
and streaming particle abundances reach about the same 
value (i.e. at ~ 80 km radius). The particle abundances 
again match well with the reference data, but the neu- 
trino spectra in the IDSA peak at lower energies than in 
the reference simulation. This is consistent with the find- 
ings above in Fig. [3] The reason is that the assumption 
of a thermal spectrum for the trapped particle compo- 
nent starts to become inaccurate at 80 km radius, while 
the trapped particle component still makes a significant 
contribution of ~ 50% to the total abundance. The in- 
accuracy of a thermal spectrum in the semi-transparent 
regim e is well-documen t ed in the supernova literature 
(e.g. iMvra &: Burrows! (|1990l )) and appears most dra- 
matic if presented on logarithmic axes. We prefer a linear 
axis where the total particle abundance, or deviations be- 
tween abundances, are proportional to the enclosed area 
between the lines. At a radius of 160 km the streaming 
particle component dominates over the trapped particle 
component as shown in panels (e) and (f). The com- 
parison of the spectra shows more accurate agreement 
because the streaming particle component keeps its full 
spectral information in the IDSA. Some deviations in the 
abundances are visible. They are most likely an effect of 
the approximations in the decoupling regime and demon- 
strate the limitations of the IDSA. 

At this stage we are curious, how much accuracy was 
lost by the decision to parameterise the trapped particle 
component by a thermal distribution function. In or- 



der to test this, we developed an alternative version that 
takes deviations from the thermal trapped particle spec- 
trum into account. The details of this test approach are 
described in Eq. (|5"B"|) in Appendix B. The resulting spec- 
tra of this more complicated version are shown in Fig. [5j 
Now the reference spectra are perfectly matched by the 
IDSA results at 40 km radius and for the antineutrinos 
also at 160 km radius. Even in the difficult regime at 
80 km radius the spectral shapes agree very nicely with 
the reference. But new deviations appear this time in 
the particle abundances, which are somewhat smaller in 
the IDSA data than in the reference data. We conclude 
from these comparisons that the results of the IDSA are 
in qualitatively robust agreement with the reference data, 
but inaccuracies of 10 — 20% are likely to be present in 
the details of the distribution function and are difficult to 
remove by simple upgrades. In order to become able to 
better distinguish the systematic errors from the acciden- 
tal errors, we include both the normal and the spectral 
version of the IDSA in the following comparisons with 
the reference data. 

Once the trapped and streaming particle components 
are specified we investigate the resulting net heating 
and cooling rates from neutrino absorption and emis- 
sion. Fig. [5] shows the energy exchange between neutri- 
nos and matter separately for the neutrinos (solid lines) 
and antineutrinos (dashed lines). This comparison has 
to be performed and interpreted very carefully, because 
the matter conditions are partially in equilibrium with 
the transported particles. Even minor deviations from 
the stationary-state equilibrium at a restart of the cal- 
culation may lead to large transient corrections in the 
exchange rates. Hence, we let the hydrodynamics, the 
reference simulation and the IDSA solution evolve after 
restart at 150 ms postbounce until the stationary-state is 
achieved in both cases. The agreement of the IDSA (lines 
labelled with 'Approx.' and 'Spectra') with the reference 
simulation (line labelled with 'Boltzmann') is excellent in 
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Fig. 4. — Particle spectra for the trapped particle component 
(dashed line), the streaming particle component (dash-dotted line), 
their sum (thin solid line) and the reference simulation (thick solid 
line) at selected radii. Panels (a) and (b) show spectra at 40 km 
radius, panels (c) and (d) at 80 km radius and panels (c) and (f) 
at 160 km radius. The panels on the left hand side (a,c,e) show 
neutrino spectra while the panels on the right hand side (b,d,f) 
show antineutrino spectra. 
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Fig. 5. — The approximate data in this figure has been calcu- 
lated with a version of the IDSA where the spectral information 
about the trapped particle component is retained as described in 
Appendix B. Otherwise, the representation is the same as in Fig. 
[4] Panels (a) and (b) show spectra at 40 km radius, panels (c) and 
(d) at 80 km radius and panels (e) and (f) at 160 km radius. The 
panels on the left hand side (a,c,e) show neutrino spectra while the 
panels on the right hand side (b,d,f) show antineutrino spectra. 



the heating region at radii larger than 120 km. Differ- 
ences of up to 30% can be seen in the cooling rates at 
smaller radii in the rates for the antineutrinos (dashed 
lines). The IDSA with a spectral representation of the 
trapped particles (lines labelled with 'Spectral') performs 
better than the IDSA with the approximate representa- 
tion of trapped particles (lines labelled with 'Approx.'). 
The reference simulation shows stronger cooling by neu- 
trinos than antineutrinos at 95 km radius. This detail 
feature is not reproduced by either IDSA version. How- 
ever, as the cooling rates are a response to an intricate 
dynamical equilibrium between the hydrodynamics and 
the quantities evolved by the transport algorithm, a sim- 
ple comparison of the heating and cooling rates alone is 
not sufficient to predict the accuracy of a long-term evo- 
lution. In analogy, the instantaneous corrections applied 
to the steering wheel of a car can vary from driver to 
driver, but most drivers follow the trajectory of the road 
on the long term. 
A comparison of profiles at the very early postbounce 



times of 1 and 3 ms is shown in Fig. [7J No significant 
deleptonisation has occured before 1 ms after bounce, 
while at 3 ms after bounce the launch of the neutrino 
burst is reflected in a trough of the electron fraction pro- 
file and a strong decline of the entropy profile at an en- 
closed mass ~ 1 Mq. Significant differences between the 
IDSA and the Boltzmann solution are only visible in the 
electron fraction profiles. 

In fact, if one defines the diffusion source as described 
in Eq. ||7J), both IDSA versions lead to an earlier and 
somewhat faster deleptonisation than in the Boltzmann 
solution. Our stationary-state assumption of the free 
streaming particle component implies that produced neu- 
trinos propagate infinitely fast. This assumption is not 
applicable in the very early postbounce phase because 
the dynamical changes of the matter conditions occur 
on a similar time scale as the neutrino transport. This 
transient deviation in the deleptonisation rates during 
the first few milliseconds lets the models start on a dif- 
ferent footing so that the performance of the approxi- 
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Fig. 6. — Comparison of the heating or cooling rates specific to 
each neutrino type at 150 ms after bounce. The solid lines show 
the energy exchange between neutrinos and matter. The dashed 
lines show the energy exchange between antineutrinos and matter. 
The lines labelled 'Approx.' and 'Spectral' have been calculated 
using the IDSA with the standard and spectral representations of 
the trapped particle component, respectively. 



mation in the more interesting later postbounce phases 
becomes difficult to assess. In order to account for the 
neutrino propagation time in an averaged way we there- 
fore limit the diffusion source E in Eq. additionally 
by fc/Ar, where Ar is an empirically determined con- 
stant parameter that stands for a typical distance the 
neutrinos have to propagate until stationary-state condi- 
tions become applicable. The value Ar = 15 km leads 
to satisfactory results. With this, the IDSA version with 
the spectral treatment of the trapped particle distribu- 
tion function (dash-dotted lines) leads to nice agreement 
with the Boltzmann results (solid lines) while the IDSA 
version with the two-parameter description of a thermal 
distribution function (dashed lines) initially deleptonises 
slightly slower, but excellent agreement is achieved later 
during the shock expansion phase. Corresponding pro- 
files at 30 ms and 100 ms after bounce are shown in Fig. 
[H The solutions based on the IDSA lead to a slightly 
more optimistic shock expansion, but the detailed fea- 
tures in the entropy and electron fraction profiles are 
accurately reproduced. 

The shock retraction phase is shown in Fig. [9l Shown 
are profiles at 150 ms and 300 ms after bounce (note 
that now the larger shock radius belongs to the earlier 
time slice). The effect of neutrino heating (and hence 
successful neutrino transport) is clearly visible in the 
entropy profile. Models that implement only neutrino 
leakage develop a much smaller entropy that is entirely 
set by the heat dissipation at the shock front. In the 
shock retraction phase, the peak entropy is obtained at 
the shock front. If neutrino heating is switched on, the 
region of peak heating is located between the neutri- 
nospheres and the shock radius (see Fig. [6]), because 
the neutrino flux dilutes at larger distances and the mat- 
ter density decreases in the outer layers. Figs. [5] and 
[5] show the negative entropy gradient between 120 km 
radius and the shock position that is caused by this neu- 
trino heating. As before we obtain excellent agreement 




Enclosed Mass [solar masses] 

Fig. 7. — Density, entropy, electron fraction and velocity as func- 
tions of enclosed mass for a model based on Boltzmann neutrino 
transport (solid lines) and two models based on the IDSA. The 
dashed lines represent data from the IDSA using a thermal trapped 
particle distribution function while the dash-dotted lines represent 
data from the IDSA with a spectral trapped particle treatment. 
The comparison is shown at two different time instances: At 1 
ms after bounce (lines with the shock discontinuity positioned be- 
tween 0.8 and 1 Mq) and at 3 ms after bounce (lines with the 
shock discontinuity positioned between 1 and 1.2 Mq). Significant 
differences are only visible in the electron fraction profiles: The 
deleptonisation in the approximate model is slightly delayed. 



in the diffusive regime. But the model with neutrino 
transport approximations starts to retract slightly ear- 
lier than the more accurate model with Boltzmann neu- 
trino transport. This is consistent with the larger cooling 
rates of the IDSA model found in Fig. [5] However, the 
overall differences are not much larger than differences 
that were also found in the first comparison of two inde- 
pend ent implementations of Boltzmann neutrino trans- 
port (|Liebend5rfer et aT]|2005f ). The chances of obtaining 
a delayed explosion seem to be more pessimistic in the 
model with approximate neutrino transport. 

A short overview of the comparison is given in Fig. 
1101 in the form of shock trajectories and neutrino lumi- 
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Fig. 8. — Same presentation of the data as in Fig. [7] Here, 
the comparison is shown at two later time instances: At 30 ms 
after bounce (lines with the shock discontinuity positioned at 150 
km radius) and at 100 ms after bounce (lines with the shock dis- 
continuity positioned at 250 km radius). The two IDSA solutions 
lead to a somewhat more optimistic shock expansion in this phase, 
but produce nicely detailed features in the electron fraction and 
entropy profiles. 



nosities for the three models. Panel (a) shows the shock 
trajectories. As stated before, the agreement is excellent 
from the early postbounce phase through the shock ex- 
pansion phase. The models with approximate neutrino 
transport initially produce a larger peak shock radius, 
but then compensate by a faster shock retraction. The 
model with the spectral treatment of the trapped parti- 
cle component (dash-dotted line) is somewhat closer to 
the reference simulation (solid line) than the standard 
IDSA with the two-parameter description of the trapped 
particle distribution function (dashed line). Panel (b) 
shows the neutrino luminosities measured in the comov- 
ing frame (this distinction is only made for the Boltz- 
mann reference results) at 500 km radius. Overall, the 
luminosities are in good agreement. The neutrino lumi- 
nosities in the approximate models are somewhat smaller 
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Fig. 9. — Same presentation of the data as in Fig. [7] Here, the 
comparison is shown at two later time instances: At 150 ms after 
bounce (lines with the shock discontinuity positioned at 230 km 
radius) and at 300 ms after bounce (lines with the shock disconti- 
nuity positioned at 150 km radius). The evolution in the diffusive 
domains is in good agreement, but the shock in the approxima- 
tive model retracts somewhat faster than in the accurate model. 
With respect to a possible delayed explosion, the approximation 
produces the more pessimistic model. 



than in the more accurate Boltzmann model. 

While experimenting with different variations in the 
approximation scheme we found that the stationary-state 
luminosity in the approximate model is not a perfectly 
smooth function of time, which leads us to the discus- 
sion of some numerical issues related to the isotropic dif- 
fusion source approach. The numerical details of our 
current implementation are outlined in Appendix B. All 
local reaction rates and the corresponding updates of the 
lepton fractions and temperature are implemented in a 
time-implicit way. The non-local contribution from the 
diffusion is also unconditionally stable because the up- 
dates arc ordered in such a way that each zone can use 
the updated distribution function of its neighbour zone in 
the upwind direction (with respect to the diffusion flux), 
while the contribution of the local distribution function 
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Fig. 10. — Panel (a) shows the shock position as a function of time for a model based on Boltzmann neutrino transport (solid line) 
and the two models using the IDSA. Both approximate models agree reasonably well with the reference model. The model based on a 
spectral treatment of the trapped particle distribution function (dash-dotted line) is somewhat closer than the standard IDSA approach 
(dashed line), but shows the same qualitative deviations of a too optimistic initial expansion and a faster shock retraction afterwards. Panel 
(b) shows the neutrino luminosities for the same three models. The agreement is satisfactory in general, the neutrino luminosities (thin 
solid lines) are generally smaller in the approximate model than in the Boltzmann model (thick solid line). The dashed lines refer to the 
antineutrino luminosities. Some experienced numerical difficulties with the approximate neutrino luminosites are discussed in the text. 



to the diffusion term is included in the time-implicit up- 
date. However, the coupling to the stationary-state so- 
lution of the streaming particle component is operator 
split. It is therefore still necessary to restrict the time 
step to obtain a numerically stable evolution. In the dif- 
fusion limit, the operator splitting is numerically stable 
due to above-mentioned implicit finite differencing. In 
the free streaming limit, the operator splitting is numer- 
ically stable due to the absence of relevant interactions 
between the streaming particle component and matter. 
It is again the semi-transparent transition regime where 
the operator splitting has the highest potential for nu- 
merical instabilities. We obtain converged solutions if 
the time step chosen is smaller than 1 km/c. This corre- 
sponds to the Courant-Friedrich-Levy (CFL) limit for the 
speed of light or the time step limit of an explicitly finite 
differenced diffusion scheme where the mean free path is 
half the zone width, i.e. in the semi-transparent regime. 
If the time step chosen is larger, the stationary-state solu- 
tion for the streaming particle component starts to jump 
from time step to time step. Rare luminosity jumps are 
still visible in Fig. [10b . All other variables, e.g. the lep- 
ton fractions and the temperature, do not react to these 
instantaneous luminosity jumps because they evolve on a 
slower time scale where the luminosity oscillations enter 
in a time-averaged manner. 

In summary, we believe that the separation of particles 
into trapped and streaming components provides an in- 
teresting ansatz for specifically tailored simplifications of 
the Boltzmann transport equations. The results in Figs. 
[Ti l 101 show that such a simplified descriptijon can repro- 
duce the key features of neutrino transport in supernova 
models. In the diffusion limit, the results are stable and 
accurate. In the transition regime the time step of our 
current implementation is at present limited by numeri- 
cal fluctuations of the stationary-state streaming particle 



flux at the transition to transparent conditions. With our 
not thoroughly optimized code, a single time step with 
the IDSA is of order 100 times faster than a correspond- 
ing time step with Boltzmann neutrino transport. As 
the fully implicit Boltzmann transport can take about 
10 times larger time steps, a full simulation with the 
IDSA is about an order of magnitude faster than the one 
with Boltzmann transport. However, this number may 
vary from application to application because the perfor- 
mance of the time-implicit Boltzmann solution does not 
scale favourably with an increased dimensionality and 
size of the computational domain. In three dimensions, 
we believe that simulations based on the IDSA are feasi- 
ble on average size high-performance computer clusters, 
while simulations with comprehensive Boltzmann trans- 
port have not yet been reported to be feasible even on 
top-performing computer systems. A possible extension 
of the IDSA to multi-dimensional applications is outlined 
in the following section. 

4. GENERALIZATION FOR MULTI-DIMENSIONAL 
APPLICATIONS 

So far, we have discussed the IDSA only in spherical 
symmetry. As mentioned above, it is not the goal of 
this paper to develop yet another approach that works 
only for spherically symmetric models. Spherical symme- 
try is much better treated with three-flavour Boltzmann 
transport where the comprehensive eq uations of radia- 
tive t r ansfer can consistently be solved (iRampp fc Jankal 
20021: iThompson et all 120031 iLiebendorfer et all l2004f 
Sumivoshi et al.ll2005f ). Hence, in this section we discuss 
how the scheme extends to the three-dimensional case. 

The state vector U of a 3D simulation only differs 
from Eq. (|17j) by the velocity v, which becomes a vec- 
tor v = (v x ,v y ,v z ). With the corresponding three direc- 
tional components i = 1 ... 3 in the momentum equation, 
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the state vector U and flux vector F become: 



U = 



V 



P 

p (e + \v 2 
pY e 
pY? 

(pziv 



I 



F = 



up 



vp 
vpvi + p 

vpY c 



vpY{ 

\ v{pz\y 



) 

(26) 

This state vector can be evolved by a standard hydrody- 
namics scheme that solves the multi-dimensional conser- 
vation law 

^-U + V ■ F = 0. (27) 

A more difficult part in the multi-dimensional IDSA is 
the consistent solution of Eqs. (J7J) and ([Tg ]) -([2"T ]) , The 
difficulty arises from the non-local scalar 



a 



-1 



3 (j + x + ■ 



(28) 



in Eq. (7}. The other equations are local and do not 
depend on the dimensionality of the problem, except that 
the integration over the angle cosine, 1/2 J dp, must be 
replaced by an integration over the entire solid angle of 
a sphere, 1/ (An) J dfl. 

We believe that there are three options when imple- 
menting Eq. (|2"5j) : (i) One may calculate Eq. (|2"5j) based 
on /* of the previous time step and update /* locally with 
Eqs. (|18 | - ([2"Tj) . Hence, only the local reactions would be 
implemented implicitly while the diffusion part remains 
explicit. This approach is only conditionally stable and 
the time step must be restricted to Ax 2 / (2 Ac). The rel- 
evance of trapped particles diminishes in the transparent 
regime around A ~ Ax. Thus, the time step restriction 
is of similar order of magnitude as the CFL condition for 
the particle propagation speed c. As the sound speed in 
neutron stars is close to the speed of light, the explicit 
approach may still be a viable option for our applica- 
tion, (ii) In order to obtain an unconditionally stable 
time step, a global implicit solution of above equations 
could also be attempted. Because of the source limit- 
ing in Eq. ([7]), the problem is most easily expressed in 
terms of the unlimited diffusion source S as function of 
df/dt as outlined in Appendix B in Eqs. |g5J) and (|33|). 
Then, one may search for a globally consistent solution 
of the field of diffusion sources. This approach has the 
disadvantage that a spectral problem has to be solved 
globally, which results in very memory- and cpu-intensive 
code. The computational effort would become compara- 
ble with multi-group flux limited diffusion or even more 
sophisticated algorithms of radiative transfer. However, 
the solution could still be restricted to the regimes where 
the trapped particle component is dominant, (iii) Per- 
haps the most promising approach is the combination of 
specific directional sweeps with a locally implicit update 
of the diffusion source, be it in th e form of a Crank- 
Nicholson scheme with ADI (e.g. iPress et al.1 (|1992l) 
or using SauPyev's asymm etric updates (|Saurvevlll96 
iTavakoli fc Davarnll I2007T ) . This approach leads to an 
unconditionally stable time step while retaining the ef- 
ficiency of an explicit scheme. In Appendix B we apply 
such an approch for our spherically symmetric applica- 
tion. 



Independent of the method used, the result of this step 
determines a partial update of the compositional quanti- 
ties in U and a spectral diffusion source S; (E) for each 
particle species / at each grid point. This information is 
stored for the following updates. The stationary-state so- 
lution for the streaming particle component is then based 



An 



- (j+x f + s 



dCt, 



(29) 



which is the natural extension of Eq. (fTU|) to the multi- 
dimensional case. The integration over d£l is again per- 
formed over the solid angle of a sphere. If the relative 
velocities in the transparent regime are large, the source 
should be transformed from the comoving frame to the 
laboratory frame (quantities measured in the laboratory 
frame carry a hat). Note that Eq. (|29|) must be solved for 
each energy group and particle species separately. But all 
directional information of the streaming particle flux VV> 
can be retrieved from the corresponding scalar potential 
ipl(E). Once tpi(E) has been calculated, the equally large 
field E; (E) is no longer needed. 
After the hydrodynamic update according to Eq. (f2"Tf , 

we have to convert the particle flux, Vipi(E), into a 
streaming particle density. This is now more ambiguous 
than in Eq. (fT3| . One may still follow the same proce- 
dure and first determine the particle scattering spheres. 
Then, one assumes a flux factor of one half in regions en- 
closed by the scattering spheres while using an isotropic 
emission ansatz along the scattering spheres to estimate 
the flux factor in the domain outside the scattering 
spheres based on their geometry. For nearly spherical 
problems such as ours, the much simpler spherically sym- 
metric ansatz described above might even be sufficient to 
obtain an estimate of the flux factor. Once the spectral 
streaming particle density 1 / (An) J ff (E) dfl has been 

determined by the quotient of Vipi (E) and the flux fac- 
tor, one cycle of updates in a multidimensional applica- 
tion is concluded. 

5. CONCLUSION 

The isotropic diffusion source approximation (IDSA) 
has been designed to treat radiative transfer most effi- 
ciently in astrophysical applications where particles dif- 
fuse out of an opaque domain that is subject to multi- 
dimensional dynamics. This scenario occurs frequently in 
new-born gravitationally bound objects, accretion discs, 
or dynamical atmospheres. We have developed and 
tested the IDSA in the context of spherically symmet- 
ric supernova models. The neutrino transport in core 
collapse supernovae provides a very challenging appli- 
cation where comprehensive solutions of the Boltzmann 
transport equation have been studied in great detail (e.g. 
iMezzacappal (|2005l ) and references therein) . Our approx- 
imation treats only the most important features of ra- 
diative transfer. In a supernova model, these are (i) the 
thermodynamics of trapped neutrinos (for example the 
dE v = —p u dV-term), (ii) the correct diffusion limit, (iii) 
the spectrum of transported neutrinos, and (iv) the angu- 
lar focussing of the neutrino propagation directions with 
increasing distance from the neutrinospheres. 

We i mplemented the IPSA in the hydrodynamics code 
Agile ijLie bcndorfe r et all [2002) to compare its perfor- 
mance with a more sophisticated transport model. We 
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compare the postbounce ev olution of a 13 M Q star 
(jNomoto fe Hashimotolfl988l ) in Newtonian gravity with 
the corresponding model N13 that is based on Boltz- 
mann neutrino transport (iLiebendorfer et alJl2005h . We 
find good agreement between the results of the IDSA 
and the reference model during the early postbounce 
phase when the neutrino burst is launched and the accre- 
tion shock expands to its maximum radius. The IDSA 
leads to a somewhat larger maximum shock radius and 
a faster shock retraction thereafter, i.e. a slightly more 
pessimistic model with respect to the possibility of a de- 
layed supernova explosion. We conclude that the concept 
of the IDSA accurately captures the feedback of neutrino 
transport onto the hydrodynamics evolution in spheri- 
cally symmetric supernova models. The neutrino fluxes 
and spectra are in nice agreement if small time steps are 
taken, but unphysical fluctuations in the stationary-state 
neutrino luminosities develop when the time steps are 
increased much beyond the neutrino propagation time 
scale. We will continue to improve the numerical stabil- 
ity of the scheme in future work. 

With respect to the choice of physical assumptions and 
with respect to the quality of the results, the IDSA is 
most closely related to the well-known flux-limited diffu- 
sion (FLD) approximation. Like the latter, the IDSA 
is guided by the diffusion limit at high opacities and 
the free streaming limit at low opacities. The FLD ap- 
proximation is based on diffusive fluxes that are limited 
by the physical particle speed at the transition to the 
free streaming regime. This limiting is well-defined in 
a spherically symmetric setup where the limiting can 
be applied to the radial velocity direction. In a multi- 
dimensional simulation, however, FLD fails to determine 
the correct direction of the net flux outside the diffu- 
sive regime. In the free streaming regime there is no 
physical connection between the local radiation intensity 
and the flux direction. A related point is the efficiency 
of the FLD approach. Since the implemented numeri- 
cal solution strategies are optimised for the diffusive do- 
main, they easily become inefficient in the free streaming 
regime, which represents often the larger part of the com- 
putational domain. These are the two main shortcom- 
ings of FLD which we address with the IDSA. The IDSA 
is based on an overlapping two-component description. 
The trapped particle component dominates in the dif- 
fusive regime while the streaming component dominates 
in the free streaming regime. In the diffusive regime, 
the IDSA is equivalent to FLD. In the semi-transparent 
regime, the IDSA is based on a different interpolation be- 
tween the diffusive and the free streaming regimes. In the 
IDSA, the limiter is applied to the sources instead of the 
fluxes. In the free streaming limit of a multi-dimensional 
application, the direction of fluxes is well-defined by the 
non-local position of the sources. The distinction be- 
tween the trapped and streaming particle components 
allows us to separately test more efficient solution algo- 



rithms for both regimes. 

For example, we assumed that the trapped particle 
component is close to thermal and weak equilibrium 
while the streaming particle component carries a non- 
equilibrium spectrum. In the framework of the IDSA we 
have compare a grey and a spectral (i.e. energy-averaged 
and energy-dependent) treatment of the evolution of the 
trapped particle component. We found visible differences 
of several percent. We believe that the differences are 
large enough to value the spectral approach, but at the 
same time they are small enough to be ignored if the ef- 
ficiency of the algorithm in 3D simulations becomes crit- 
ical. However, we emphasise that this comparison must 
not be confused with the comparison of a globally grey 
scheme and a globally spectral approach. We believe that 
the streaming particle component, which represents the 
lepton and energy flux in a supernova application, must 
be treated in a spectral manner. The location of the 
neutrinospheres is highly energy-dependent. Hence, the 
geometry of the radiative transfer problem is different for 
each neutrino energy. An energy-averaged local quantity 
can only be expected to be accurately represented if one 
averages over solutions of the energy-dependent trans- 
port problem. A grey scheme, that solves an energy- 
averaged transport problem instead, is unlikely to pro- 
duce accurate results in the supernova context. 

Aside from the field of neutrino-driven supernova 
models, FLD has been used in fields ranging f rom 
planet migr ation ( Paardckoopcr fc Mellcnia |20(M ) to 
cosmology (jBonanno fe Romano! Il994h . includin g ar- 
eas s uch as accretion discs, both around stars (jKlevi 
1989) and more energetic objects where the accretion 
disc is dominated by radiation (e.g. (jTurner et al.l 
l2002t iHuieirat fe Camenzindll200Clf )). and star formation 
(jWhitehouse fc Batd l2006) . The IDSA is potentially ap- 
plicable to these and other areas of astrophysics where 
FLD is commonly used. However, neither the IDSA nor 
the FLDA can replace the development of more reliable 
transport schemes because, for each application, it has 
to be ascertained that no important physical ingredients 
are missed. We hope that the IDSA is well-suited to ac- 
company the application of sophisticated radiative trans- 
fer algorithms in order to enable an efficient exploration 
of input physics, parameter space, or dimensionality in 
computer models where a comprehensive solution of the 
Boltzmann transport equation is not affordable. 
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APPENDIX A: THE DIFFUSION LIMIT REVISITED 

Even if the following is old wine in a new hose, we derive the diffusion limit here in brief in order to make a clear 
distinction between a flux-limited diffusion approximation, which determines the particle flux based on the diffusion 
limit and our scheme, which requires a well-founded approximation for the collision integral in the diffusive limit. 

The source term on the right hand side of the Boltzmann equation ([3]) features an isotropic emissivity, j, an isotropic 
absorptivity, and an isoenergetic scattering kernel, R. For the derivation of the diffusion l i mit, the angular depen- 
dence of the scattering kernel is usually expanded in a Legendre series (see e.g. iBruennl ()1985l ): iMezzacappa &: Bruennl 
(Gil)), 

-^-3 R (E, ft, = T" E (21 + 1) <t>t (E) f Pi (cos 9) dip, (30) 
c (he) 47r £ Jo 

1 /2 

where cos 8 = /i/i' + cos ip [(l — /i 2 ) (l — m' 2 )] specifies the angle between the incoming and the scattered particle 
propagation directions. In Eq. (|30"|) we define (j>i (E) such that it enters Eq. ([3]) without additional factors. If we 
denote the operator on the left hand side of Eq. (J3]) by D (/) and truncate the Legendre expansion after the second 
term, we reach a very concise representation of the Boltzmann equation, 

D (/) =]-(]+X + Mf + 4>v\ J + fdfi + Zft<t>i\ J + ffidfx. (31) 

From the zeroth and first angular moments of Eq. (|31|) we extract an expression for the zeroth and first moments of 
the distribution function, respectively, 
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If we substitute these terms back into Eq. 
1 



/ 



3 + X + 0o 



3 -D(f) + 



we can express the distribution function / in the following way, 



3 + X 



3 



3 + X + <t>o - <t> 



-i f D(f)vdfi 



(33) 



Now, we perform a Chapman-Enskog expansion (Chapm an fc CowlingllT97Clf L To this purpose we introduce a small 
expansion parameter, e, and replace the large source functions j, \ an d 0o,i hi Eq. ()33j) by j/e, x/ £ an d < t'o,i/^i 
respectively. From the zeroth order terms in e we then obtain as expected the isotropic equilibrium distribution 
function fo — j / (j + \)- However, we can also determine the first order corrections to the equilibrium distribution 
function: 



e.h = 



-1 



3 + X + 4>o 



D(fo) 



3+X% J 



3 + X + ^0 



(34) 



After these preparations, we obtain the leading order particle number exchange rate with matter, s, by substituting 
the approximate distribution function / ~ / + ef\ into Eq. (|3Tj) and integrating the left hand side over particle 
propagation angles: 

s= l -fD{f + eh)d l i. (35) 

This integral can be simplified if one decomposes the operator D (f) into an operator D + (f) containing terms that are 
symmetric in \i and an operator D~ (f) containing terms that are antisymmetric in \i. Hence, D (/) = D + (f)+D~ (f), 
with 
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+ ff\ d f , ( dln P , 3w\ H 2 \ 9 f , [ 2 fdlnp 3v\ 
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dE 



- (1 - M 2 ) — 



(36) 



Only the symmetric terms D + (fo), D + [D + (fo)] and D~ [D~ (fo)] survive the angular integration in Eq. 
Furthermore, we may neglect the D + [D + (/)] terms because they are of higher order in the expansion than the 
leading 1/2 J D + (fo) d[i term. The D~ [D~ (fo)] terms, however, are of leading order because there is no lower order 
contribution involving D~ . The intermediate result 
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can be further simplified if one explicitly calculates the integral in the last term using D (fo) 
us with the following leading order source terms in the diffusion limit, 



(fo) d ^-\ 
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-D- (fo) df, 



fidfo/dr. This leaves 



(38) 



The original operators in Eq. (|36[) may now be readily substituted. Several c'/ctyi-terms drop out due to the isotropy 
of fo = j I (j + x) an d the remaining angular integrations can easily be performed. With this, we obtain the following 
simple expression for the collision integral in the diffusion limit, 



dfo 
cdt 



1 dlnp E dfo 
3 cdt dE 



J_d_ 

r 2 dr 



dfo 



3 (j + X + 4>) dr 



(39) 



where, in order to shorten the notation, we define = <f>o — 4>i ■ Note that a flux-limited diffusion scheme applies Eq. 
(|39[) as evolution equation for an unknown / (t, r, E) (instead of fo) together with a source function s (f). Our scheme 
uses Eq. (f3"9"]) as leading order approximation to the source s for a given equilibrium distribution function fo, exactly 
as derived. The first term in Eq. ([39]) describes sources that arise from adjustments to a time-dependent equilibrium, 
the second term accounts for changes in the average particle energy when the fluid is compressed or expanded and the 
third term describes the divergence of the diffusive flux in the rest frame of the fluid according to the transport mean 
free path X t = (j + x + (t>)~~ ■ 

APPENDIX B: IMPLEMENTATION AND FINITE DIFFERENCING 

In this appendix we discuss the finite differencing of our implementation of the IDSA. Our scheme solves Eq. (|18[) 
for the trapped particles, Eq. (jlOp for the streaming particles, and Eq. (|14p for the hydrodynamics in an operator 
splitting manner. The spherically symmetric computational domain is represented by an array of concentric fluid 
shells labelled by the index i — 1 . . . i rnax . Unprimed indices are used for zone centre values and primed indices for the 
upper zone edge values. Energy-dependent quantities are calculated for a set of discrete particle energies Ek, where 
k = 1 . . . k max . This representation exists for several time instances, which we label by an upper index n = 1 . . . n max . 
For example, A" k refers to the mean free path A (r i+1 / 2 , Ek, t n ) at the zone edge between r% and fi+i, at energy Ek 
and at time t n . 
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A time steps begins with the construction of the particle flux at time t n+1 from the particle sources, X — 
(j + X) h J f s d[J-, which have been calculated and stored during the previous time step t n . The particle flux is given by 
the gradient of the potential ip as described in Eq. (|10p . In spherical symmetry, this Poisson equation is easily solved 
by the integration of the sources over the volume, 



dip 
dr 



n+l 
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fa 



(40) 



where we used a tilde to denote stimulated absorption x = j + X- At present we neglect the Lorentz transformation 
of the sources from the comoving frame to the laboratory frame. 

Next, we retrieve the current conditions from the solution vector U" = U(ri,t n ) and calculate the energy-dependent 
transport mean free path A" fc and the optical depth r" fc = (1/A(r)^) dr. In spherical symmetry, the energy- 
dependent radius R™ k of the neutrinosphere is found at the point where T^i?" fc ^ = 2/3. The neutrinosphere is 

then used to derive the streaming particle density from the streaming particle flux according to Eq. (|13p . Again, 
we neglect the Lorentz transformation between the laboratory and comoving frames and approximate the streaming 
particle density by 
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(41) 



The factor (ry-i/rj) converts the flux from the inner zone edge %' — 1, where it has been calculated in Eq. (|40|) . to 
the zone centre 7"j, where it is used in the calculation of the diffusion source. Two measures significantly increase the 
stability of the scheme. On the one hand it is important to exclude the contribution of zone i to the streaming particle 
flux. On the other hand, the maximum function with the square brackets eliminates fluxes that stream against the 
density gradient. 

The trapped particle distribution functions are separated into a thermal component, and a spectral perturbation, 
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The thermal part is represented by a Fermi function, /^ h , where the particle temperature and particle degeneracy are 
chosen such that Eq. P^|) is fulfilled for the current values of the quantities Y i ' n and Z\ n in the solution vector U™. If 
the analytical approximations given in (|Epstein fc Pethicklll98ll) are used, a consistent set of the particle temperature 
and degeneracy parameters can be obtained from the solution of a one-dimensional Newton-Raphson scheme. For the 
standard version of the IDSA we neglect the perturbations by setting Sfj'u = 0. In order to test this approximation 
in Sect. [31 we implemented a second version of the IDSA, where the deviation of the trapped particle spectrum from 
the thermal spectrum is evolved in time as described in Eq. (|55|) . 

The next step is the search for a consistent solution of Eqs. (fT5)) . (f2"0"| and (|2"Tj) . In most of the following monochro- 
matic equations, all quantities carry the same energy index k. To further reduce the amount of indices we will write 
the energy index only if quantities at different energies are involved in one equation. We discretise the fast reaction 
rates on the right hand side of Eq. (|18[) in a time-implicit way, 
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and eliminate the dependence on f i on the right hand side by rewriting Eq. (|43|) in the form 
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The diffusion source S, which is defined in Eq. ([?])■ depends on a. We will work with the following definition of a 
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(48) 



In the search for a consistent diffusion source we first assume that S = a + x\ J / S ^M i s n °t limited. Then we 
apply the limiters in Eq. J7|) to the result. It is known that the numerical evolution of the diffusion equation is not 
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uncon ditionally stable unless the term a is finite differenced in an at least partially implicit manner (e.g. I Press et al.l 
(1992)). Hence, we would like to discretise the diffusion source as 
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(49) 



The dependence of this equation on the unknown distribution functions f t ' n+1 can be removed if we substitute Eq. 
(|4"4"|) in the first three terms. The result is a non-local equation involving S£i , and , which could globally 

be solved for However, since our target application is based on a three-dimensional parallelised hydrodynamics 

code, for the moment we try to avoid a global solution of Eqs. (149)) and (f44|) . At least in our spherically symmetric 
example, the diffusive fluxes propagate almost exclusively outward. Hence we finite difference the outward propagating 
flux implicitly and the inward propagating flux explicitly, which leads to 
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If one applies Eq. (|50|) from the inside out, is known from the update of the previous zone. After substitution 

of Eq. (|4"4")) in the first term, Eq. (|50")) can be solved for so that the limiter from Eq. ([7]) can be applied, 
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(51) 



This diffusion source can now be used in Eq. (|4"4"|) to determine the discrete time derivative of the trapped 

particle distribution function, f/^'™^ 1 — / (cA<). Furthermore, this derivative specifies the net interaction rate 
between matter and trapped particles, which can be calculated according to Eq. 



f 



t,n+l 



f 



cAt 



/■s.n+1 i 



(52) 



In the application of Eq. (|5"Tj) to the supernova model in Sect. [3] we observed a slightly too rapid deleptonisation 
few milliseconds after the bounce of the stellar core. As the dynamical time scale of the bounce is comparable to the 
neutrino propagation time scale in this short transition phase, the stationary-state assumption of our approximation 
may not hold. This effect disappeared when we additionally limited S in Eq. (|5ip by / t,n+1 /Ar, where Ar ~ 15 km 
is an empirically determined constant parameter. 

In the case of neutrino transport with electron neutrinos and electron antineutrinos Eq. (|52p leads to the following 
updates of the electron fraction and the internal specific energy in Eqs. (f2"0"l) and ([2T]) : 



Y 



n+l 



yr. 



cAt 

n+l 



cAt 



mb 4-7T 
mb 47r 

3 a 



S^Sfe = 



P? (he) 



and x(p?,K 



n+l 



(53) 
(54) 

which enter the right 



However, the neutrino-matter interactions j" 

hand sides of Eqs. ([53)1 and ((54)) via Eq. (|44| . depend on the updated values for the electron fraction and specific 
internal energy. Hence we solve this nonlinear algebraic system of equations numerically by two-dimensional Newton- 
Raphson iterations to determine a consistent j^™" 1 " 1 , e™ +1 )-pair 1 . This consistency is important to accurately represent 
reactive equilibria at high opacities. Following an approach of Mczzacap pa &: Messerl (|1999l ). we first evaluate Eq. (|52|) 
on four corners of a rectangle, centred around the initial values {^.™,e™} in the space spanned by electron fraction 
and specific energy. Then, we use interpolations in this rectangle to efficiently obtain the reaction rates and their 
derivatives with respect to Y c and e during the iterations required to find a consistent solution of Eqs. (|53|) and (|54| . 
Once the iterations have converged, we update pY e and pe in the hydrodynamics state vector U n , where n! indicates 
an intermediate incomplete update. 

Finally, the components pY t and (pZ 4 ) 3 ^ 4 in U n are updated based on Eqs. (|22"|) . ([23)1 and ((44)) . In the standard 
version of the IDSA, the spectral information of the trapped particle component is discarded at this point. For the 



In fact, for compatibility reasons with the hydrodynamics code and the equation of state, our implementation of the Newton-Raphson 
scheme works with the logarithmic entropy instead of the specific energy. 
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alternative version, the spectral information is stored in an energy-dependent vector 5fl' k +l that will be used in the 
next time step according to Eq. (|^|) . The latter is obtained from 

5 fk=fi,k -Ik (Ji , Z i ) 

Xf t.n+i_cf* , a oh - a±bi - (a bi - aib ) E k th / t .n+i ~t,ra+A , KK s 
A/i,* (6? _ ^ /* ^ ,^ J (55) 



b m = 



£/f (Y^ + \Z^ +1 )E^dE k . 



The first line in Eq. (|55p calculates the deviation between the updated trapped particle distribution function and the 
thermal spectrum. Due to numerical inaccuracies, spurious errors can accumulate in 6f£ over many time steps and start 
to contribute to the trapped particle abundance and the mean specific energy Z\ which should be determined by the 

thermal part alone. Adding the correction term on the second line of Eq. ([55]) guarantees that <5 fi'k +1 E k dE k = 
and J2 &flt +lE l dE k = to machine precision. The correction term depends on energy moments of 5f£ and the 
thermal distribution function. 

The cycle of updates is concluded by evolving the partially updated hydrodynamics state U n to U n+1 according to 
Eq. (I14p (and including gravity) by any standard hydro dynamics scheme. For th is demonstration of the method we 
use the spherically symmetric hydrodynamics code Agile (|Liebendorfer et aLll2002f ) with a second order accurate TVD 
advection scheme. 



